!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine weno (ro, mx, my, mz, bx, by, bz, en, az, po, gx, gy, x, y,                                                             &
   nx, ny, ng, n_neighs, dx, dy, dt, max_et, gm, cd, cfl, div_clean_flag, resistivity_flag, riemann_solver_flag,                   &
   fxro_bnd, fxmx_bnd, fxmy_bnd, fxmz_bnd, fxbx_bnd, fxby_bnd, fxbz_bnd, fxen_bnd, fxpo_bnd,                                       &
   fyro_bnd, fymx_bnd, fymy_bnd, fymz_bnd, fybx_bnd, fyby_bnd, fybz_bnd, fyen_bnd, fypo_bnd)
!==================================================================================================================================|

   use interface_check_state
   implicit none

   integer(4), intent(in) :: nx, ny, ng, n_neighs, div_clean_flag, resistivity_flag, riemann_solver_flag
   real(8), intent(in) :: dx, dy, dt, gm, cd, cfl, x(nx), y(ny)
   real(8), intent(inout) :: max_et
   real(8), dimension(nx, ny), intent(inout) :: ro, mx, my, mz, bx, by, bz, en, az, po, gx, gy
   real(8), dimension(ny, n_neighs), intent(inout) ::                                                                              &
      fxro_bnd, fxmx_bnd, fxmy_bnd, fxmz_bnd, fxbx_bnd, fxby_bnd, fxbz_bnd, fxen_bnd, fxpo_bnd
   real(8), dimension(nx, n_neighs), intent(inout) ::                                                                              &
      fyro_bnd, fymx_bnd, fymy_bnd, fymz_bnd, fybx_bnd, fyby_bnd, fybz_bnd, fyen_bnd, fypo_bnd

   integer(4) :: i, im1, j, jm1
   real(8), dimension(nx, ny) :: ro_pre, mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, en_pre, az_pre, po_pre
   real(8), dimension(nx, ny) :: fxro, fxmx, fxmy, fxmz, fxbx, fxby, fxbz, fxen, fxpo
   real(8), dimension(nx, ny) :: fyro, fymx, fymy, fymz, fybx, fyby, fybz, fyen, fypo
   real(8), dimension(nx, ny) :: bx_resis, by_resis, bz_resis, en_resis, az_resis !, sx, sy, te_cond_x, te_cond_y

   real(8) :: ch, ch2, cp2
   real(8) :: eps, gmm1  !, vx, vy, vz, b2, v2, pr, c2, s2, ca2, cfx, cfy
   real(8) :: dt_dx, dt_dy, dx2, dy2

!----------------------------------------------------------------------------------------------------------------------------------|
   gmm1 = gm - 1.0d0
   eps = 1.0d-12
   dx2 = dx * 2.0d0
   dy2 = dy * 2.0d0
   dt_dx = dt / dx
   dt_dy = dt / dy
   ch = cfl * min (dx, dy) / dt * div_clean_flag
   ch2 = ch * ch
!   cp2 = ch * cd / 3.0d0 + eps
   cp2 = -dt * ch2 / log (cd) + eps

!----------------------------------------------------------------------------------------------------------------------------------|
!  first step of 2nd order Runge Kutta method
!----------------------------------------------------------------------------------------------------------------------------------|
!   do j = 1, ny
!   do i = 1, nx
!      vx = mx(i, j) / ro(i, j) 
!      vy = my(i, j) / ro(i, j) 
!      vz = mz(i, j) / ro(i, j) 
!      b2 = bx(i, j) * bx(i, j) + by(i, j) * by(i, j) + bz(i, j) * bz(i, j)
!      v2 = vx * vx + vy * vy + vz * vz
!      pr = (en(i, j) - v2 * ro(i, j) / 2.0d0 - b2 / 2.0d0) * gmm1
!      c2 = gm * pr
!      s2 = c2 + b2
!      ca2 = bx(i, j) * bx(i, j)
!      cfx = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * ca2)) / ro(i, j)  / 2.0d0)
!      sx(i, j) = cfx + abs (vx)
!      ca2 = by(i, j) * by(i, j)
!      cfy = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * ca2)) / ro(i, j)  / 2.0d0)
!      sy(i, j) = cfy + abs (vy)
!   enddo
!   enddo

   call weno_lx (fxro, fxmx, fxmy, fxmz, fxbx, fxby, fxbz, fxen, fxpo,                                                             &
      ro, mx, my, mz, bx, by, bz, en, po, nx, ny, gm, ch2, riemann_solver_flag)
   call weno_ly (fyro, fymx, fymy, fymz, fybx, fyby, fybz, fyen, fypo,                                                             &
      ro, mx, my, mz, bx, by, bz, en, po, nx, ny, gm, ch2, riemann_solver_flag)

   do j = 2, ny
   do i = 2, nx
      im1 = i - 1
      jm1 = j - 1
      ro_pre(i, j) = ro(i, j) - dt_dx * (fxro(i, j) - fxro(im1, j)) - dt_dy * (fyro(i, j) - fyro(i, jm1))
      mx_pre(i, j) = mx(i, j) - dt_dx * (fxmx(i, j) - fxmx(im1, j)) - dt_dy * (fymx(i, j) - fymx(i, jm1))
      my_pre(i, j) = my(i, j) - dt_dx * (fxmy(i, j) - fxmy(im1, j)) - dt_dy * (fymy(i, j) - fymy(i, jm1))
      mz_pre(i, j) = mz(i, j) - dt_dx * (fxmz(i, j) - fxmz(im1, j)) - dt_dy * (fymz(i, j) - fymz(i, jm1))
      bx_pre(i, j) = bx(i, j) - dt_dx * (fxbx(i, j) - fxbx(im1, j)) - dt_dy * (fybx(i, j) - fybx(i, jm1))
      by_pre(i, j) = by(i, j) - dt_dx * (fxby(i, j) - fxby(im1, j)) - dt_dy * (fyby(i, j) - fyby(i, jm1))
      bz_pre(i, j) = bz(i, j) - dt_dx * (fxbz(i, j) - fxbz(im1, j)) - dt_dy * (fybz(i, j) - fybz(i, jm1))
      en_pre(i, j) = en(i, j) - dt_dx * (fxen(i, j) - fxen(im1, j)) - dt_dy * (fyen(i, j) - fyen(i, jm1))
      po_pre(i, j) = po(i, j) - dt_dx * (fxpo(i, j) - fxpo(im1, j)) - dt_dy * (fypo(i, j) - fypo(i, jm1))
      az_pre(i, j) = az(i, j)
   enddo
   enddo

   call check_state (ro_pre, mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, en_pre, x, y, nx, ny, gm)

   if (resistivity_flag .eq. 0) then
      do j = 1, ny
      do i = 1, nx
         bx_resis(i, j) = 0.0d0
         by_resis(i, j) = 0.0d0
         bz_resis(i, j) = 0.0d0
         en_resis(i, j) = 0.0d0
         az_resis(i, j) = 0.0d0
      enddo
      enddo
      max_et = 0.0d0
   else
      call resistivity (bx_resis, by_resis, bz_resis, en_resis, az_resis,                                                          &
         ro, mx, my, mz, bx, by, bz, en, x, y, dx, dy, nx, ny, max_et)
   endif

!   call thermal_conduction (te_cond_x, te_cond_y, ro, mx, my, mz, bx, by, bz, en, thermal_amp, gm, dx, dy, nx, ny)

   call source_term (mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, en_pre, az_pre, po_pre,                                       &
      ro, mx, my, bx, by, bz, po, gx, gy, bx_resis, by_resis, bz_resis, en_resis, az_resis,                                        &
      ch2, cp2, dx, dy, dt, nx, ny, div_clean_flag)

!----------------------------------------------------------------------------------------------------------------------------------|
!  bnd value in x direction
!----------------------------------------------------------------------------------------------------------------------------------|
   do j = 1, ny
      fxro_bnd(j, 1) = fxro(ng, j) * 0.5d0
      fxro_bnd(j, 2) = fxro(nx - ng, j) * 0.5d0

      fxmx_bnd(j, 1) = fxmx(ng, j) * 0.5d0
      fxmx_bnd(j, 2) = fxmx(nx - ng, j) * 0.5d0

      fxmy_bnd(j, 1) = fxmy(ng, j) * 0.5d0
      fxmy_bnd(j, 2) = fxmy(nx - ng, j) * 0.5d0

      fxmz_bnd(j, 1) = fxmz(ng, j) * 0.5d0
      fxmz_bnd(j, 2) = fxmz(nx - ng, j) * 0.5d0

      fxbx_bnd(j, 1) = fxbx(ng, j) * 0.5d0
      fxbx_bnd(j, 2) = fxbx(nx - ng, j) * 0.5d0

      fxby_bnd(j, 1) = fxby(ng, j) * 0.5d0
      fxby_bnd(j, 2) = fxby(nx - ng, j) * 0.5d0

      fxbz_bnd(j, 1) = fxbz(ng, j) * 0.5d0
      fxbz_bnd(j, 2) = fxbz(nx - ng, j) * 0.5d0

      fxen_bnd(j, 1) = fxen(ng, j) * 0.5d0
      fxen_bnd(j, 2) = fxen(nx - ng, j) * 0.5d0

      fxpo_bnd(j, 1) = fxpo(ng, j) * 0.5d0
      fxpo_bnd(j, 2) = fxpo(nx - ng, j) * 0.5d0
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  bnd value in y direction
!----------------------------------------------------------------------------------------------------------------------------------|
   do i = 1, nx
      fyro_bnd(i, 1) = fyro(i, ng) * 0.5d0
      fyro_bnd(i, 2) = fyro(i, ny - ng) * 0.5d0

      fymx_bnd(i, 1) = fymx(i, ng) * 0.5d0
      fymx_bnd(i, 2) = fymx(i, ny - ng) * 0.5d0

      fymy_bnd(i, 1) = fymy(i, ng) * 0.5d0
      fymy_bnd(i, 2) = fymy(i, ny - ng) * 0.5d0

      fymz_bnd(i, 1) = fymz(i, ng) * 0.5d0
      fymz_bnd(i, 2) = fymz(i, ny - ng) * 0.5d0

      fybx_bnd(i, 1) = fybx(i, ng) * 0.5d0
      fybx_bnd(i, 2) = fybx(i, ny - ng) * 0.5d0

      fyby_bnd(i, 1) = fyby(i, ng) * 0.5d0
      fyby_bnd(i, 2) = fyby(i, ny - ng) * 0.5d0

      fybz_bnd(i, 1) = fybz(i, ng) * 0.5d0
      fybz_bnd(i, 2) = fybz(i, ny - ng) * 0.5d0

      fyen_bnd(i, 1) = fyen(i, ng) * 0.5d0
      fyen_bnd(i, 2) = fyen(i, ny - ng) * 0.5d0

      fypo_bnd(i, 1) = fypo(i, ng) * 0.5d0
      fypo_bnd(i, 2) = fypo(i, ny - ng) * 0.5d0
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  second step of 2nd order Runge Kutta method
!----------------------------------------------------------------------------------------------------------------------------------|
!   do j = 1, ny
!   do i = 1, nx
!      vx = mx_pre(i, j) / ro_pre(i, j) 
!      vy = my_pre(i, j) / ro_pre(i, j) 
!      vz = mz_pre(i, j) / ro_pre(i, j) 
!      b2 = bx_pre(i, j) * bx_pre(i, j) + by_pre(i, j) * by_pre(i, j) + bz_pre(i, j) * bz_pre(i, j)
!      v2 = vx * vx + vy * vy + vz * vz
!      pr = (en_pre(i, j) - v2 * ro_pre(i, j) / 2.0d0 - b2 / 2.0d0) * gmm1
!      c2 = gm * pr
!      s2 = c2 + b2
!      ca2 = bx_pre(i, j) * bx_pre(i, j)
!      cfx = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * ca2)) / ro_pre(i, j)  / 2.0d0)
!      sx(i, j) = cfx + abs (vx)
!      ca2 = by_pre(i, j) * by_pre(i, j)
!      cfy = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * ca2)) / ro_pre(i, j)  / 2.0d0)
!      sy(i, j) = cfy + abs (vy)
!      if (isnan(sx(i, j)) .or. sx(i, j) .gt. 1.0d0 / eps .or. sx(i, j) .lt. -1.0d0 / eps) sx(i, j) = 0.0d0
!      if (isnan(sy(i, j)) .or. sy(i, j) .gt. 1.0d0 / eps .or. sy(i, j) .lt. -1.0d0 / eps) sy(i, j) = 0.0d0
!   enddo
!   enddo

   call weno_lx (fxro, fxmx, fxmy, fxmz, fxbx, fxby, fxbz, fxen, fxpo,                                                             &
      ro_pre, mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, en_pre, po_pre, nx, ny, gm, ch2, riemann_solver_flag)
   call weno_ly (fyro, fymx, fymy, fymz, fybx, fyby, fybz, fyen, fypo,                                                             &
      ro_pre, mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, en_pre, po_pre, nx, ny, gm, ch2, riemann_solver_flag)

   do j = 2, ny
   do i = 2, nx
      im1 = i - 1
      jm1 = j - 1
      ro(i, j) = (ro(i, j) + ro_pre(i, j) - dt_dx * (fxro(i, j) - fxro(im1, j)) - dt_dy * (fyro(i, j) - fyro(i, jm1))) * 0.5d0
      mx(i, j) = (mx(i, j) + mx_pre(i, j) - dt_dx * (fxmx(i, j) - fxmx(im1, j)) - dt_dy * (fymx(i, j) - fymx(i, jm1))) * 0.5d0
      my(i, j) = (my(i, j) + my_pre(i, j) - dt_dx * (fxmy(i, j) - fxmy(im1, j)) - dt_dy * (fymy(i, j) - fymy(i, jm1))) * 0.5d0
      mz(i, j) = (mz(i, j) + mz_pre(i, j) - dt_dx * (fxmz(i, j) - fxmz(im1, j)) - dt_dy * (fymz(i, j) - fymz(i, jm1))) * 0.5d0
      bx(i, j) = (bx(i, j) + bx_pre(i, j) - dt_dx * (fxbx(i, j) - fxbx(im1, j)) - dt_dy * (fybx(i, j) - fybx(i, jm1))) * 0.5d0
      by(i, j) = (by(i, j) + by_pre(i, j) - dt_dx * (fxby(i, j) - fxby(im1, j)) - dt_dy * (fyby(i, j) - fyby(i, jm1))) * 0.5d0
      bz(i, j) = (bz(i, j) + bz_pre(i, j) - dt_dx * (fxbz(i, j) - fxbz(im1, j)) - dt_dy * (fybz(i, j) - fybz(i, jm1))) * 0.5d0
      en(i, j) = (en(i, j) + en_pre(i, j) - dt_dx * (fxen(i, j) - fxen(im1, j)) - dt_dy * (fyen(i, j) - fyen(i, jm1))) * 0.5d0
      po(i, j) = (po(i, j) + po_pre(i, j) - dt_dx * (fxpo(i, j) - fxpo(im1, j)) - dt_dy * (fypo(i, j) - fypo(i, jm1))) * 0.5d0
   enddo
   enddo

   if (resistivity_flag .eq. 0) then
      do j = 1, ny
      do i = 1, nx
         bx_resis(i, j) = 0.0d0
         by_resis(i, j) = 0.0d0
         bz_resis(i, j) = 0.0d0
         en_resis(i, j) = 0.0d0
         az_resis(i, j) = 0.0d0
      enddo
      enddo
      max_et = 0.0d0
   else
      call resistivity (bx_resis, by_resis, bz_resis, en_resis, az_resis,                                                          &
         ro_pre, mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, en_pre, x, y, dx, dy, nx, ny, max_et)
   endif

!   call thermal_conduction (te_cond_x, te_cond_y, ro_pre, mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, en_pre,                  &
!      thermal_amp, gm, dx, dy, nx, ny)

   call source_term (mx, my, mz, bx, by, bz, en, az, po,                                                                           &
      ro_pre, mx_pre, my_pre, bx_pre, by_pre, bz_pre, po_pre, gx, gy,                                                              &
      bx_resis, by_resis, bz_resis, en_resis, az_resis,                                                                            &
      ch2, cp2, dx, dy, dt / 2.0d0, nx, ny, div_clean_flag)

!----------------------------------------------------------------------------------------------------------------------------------|
!  bnd value in x direction
!----------------------------------------------------------------------------------------------------------------------------------|
   do j = 1, ny
      fxro_bnd(j, 1) = fxro_bnd(j, 1) + fxro(ng, j) * 0.5d0
      fxro_bnd(j, 2) = fxro_bnd(j, 2) + fxro(nx - ng, j) * 0.5d0

      fxmx_bnd(j, 1) = fxmx_bnd(j, 1) + fxmx(ng, j) * 0.5d0
      fxmx_bnd(j, 2) = fxmx_bnd(j, 2) + fxmx(nx - ng, j) * 0.5d0

      fxmy_bnd(j, 1) = fxmy_bnd(j, 1) + fxmy(ng, j) * 0.5d0
      fxmy_bnd(j, 2) = fxmy_bnd(j, 2) + fxmy(nx - ng, j) * 0.5d0

      fxmz_bnd(j, 1) = fxmz_bnd(j, 1) + fxmz(ng, j) * 0.5d0
      fxmz_bnd(j, 2) = fxmz_bnd(j, 2) + fxmz(nx - ng, j) * 0.5d0

      fxbx_bnd(j, 1) = fxbx_bnd(j, 1) + fxbx(ng, j) * 0.5d0
      fxbx_bnd(j, 2) = fxbx_bnd(j, 2) + fxbx(nx - ng, j) * 0.5d0

      fxby_bnd(j, 1) = fxby_bnd(j, 1) + fxby(ng, j) * 0.5d0
      fxby_bnd(j, 2) = fxby_bnd(j, 2) + fxby(nx - ng, j) * 0.5d0

      fxbz_bnd(j, 1) = fxbz_bnd(j, 1) + fxbz(ng, j) * 0.5d0
      fxbz_bnd(j, 2) = fxbz_bnd(j, 2) + fxbz(nx - ng, j) * 0.5d0

      fxen_bnd(j, 1) = fxen_bnd(j, 1) + fxen(ng, j) * 0.5d0
      fxen_bnd(j, 2) = fxen_bnd(j, 2) + fxen(nx - ng, j) * 0.5d0

      fxpo_bnd(j, 1) = fxpo_bnd(j, 1) + fxpo(ng, j) * 0.5d0
      fxpo_bnd(j, 2) = fxpo_bnd(j, 2) + fxpo(nx - ng, j) * 0.5d0
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  bnd value in y direction
!----------------------------------------------------------------------------------------------------------------------------------|
   do i = 1, nx
      fyro_bnd(i, 1) = fyro_bnd(i, 1) + fyro(i, ng) * 0.5d0
      fyro_bnd(i, 2) = fyro_bnd(i, 2) + fyro(i, ny - ng) * 0.5d0

      fymx_bnd(i, 1) = fymx_bnd(i, 1) + fymx(i, ng) * 0.5d0
      fymx_bnd(i, 2) = fymx_bnd(i, 2) + fymx(i, ny - ng) * 0.5d0

      fymy_bnd(i, 1) = fymy_bnd(i, 1) + fymy(i, ng) * 0.5d0
      fymy_bnd(i, 2) = fymy_bnd(i, 2) + fymy(i, ny - ng) * 0.5d0

      fymz_bnd(i, 1) = fymz_bnd(i, 1) + fymz(i, ng) * 0.5d0
      fymz_bnd(i, 2) = fymz_bnd(i, 2) + fymz(i, ny - ng) * 0.5d0

      fybx_bnd(i, 1) = fybx_bnd(i, 1) + fybx(i, ng) * 0.5d0
      fybx_bnd(i, 2) = fybx_bnd(i, 2) + fybx(i, ny - ng) * 0.5d0

      fyby_bnd(i, 1) = fyby_bnd(i, 1) + fyby(i, ng) * 0.5d0
      fyby_bnd(i, 2) = fyby_bnd(i, 2) + fyby(i, ny - ng) * 0.5d0

      fybz_bnd(i, 1) = fybz_bnd(i, 1) + fybz(i, ng) * 0.5d0
      fybz_bnd(i, 2) = fybz_bnd(i, 2) + fybz(i, ny - ng) * 0.5d0

      fyen_bnd(i, 1) = fyen_bnd(i, 1) + fyen(i, ng) * 0.5d0
      fyen_bnd(i, 2) = fyen_bnd(i, 2) + fyen(i, ny - ng) * 0.5d0

      fypo_bnd(i, 1) = fypo_bnd(i, 1) + fypo(i, ng) * 0.5d0
      fypo_bnd(i, 2) = fypo_bnd(i, 2) + fypo(i, ny - ng) * 0.5d0
   enddo

!----------------------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine weno
